##########################################################
### Estimate First Stage for Control Function Approach ###
##########################################################
	
cluster_run <- FALSE # Run on UW cluster	
network_flag <- FALSE
churn_flag <- FALSE	
	
##### Specify covariates

	add_year_FEs <- TRUE
	add_market_FEs <- TRUE
	
	product_covariates <- c("Premium","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net","Anthem_HMO")
	instrumental_variables <- c("geog_factor","Hausman_instrument","HMO_own","HMO_other","AV_own","AV_other")
	
##### Part 1: Perform First Stage at the Plan Level
	
# Load Data 

	if(cluster_run) {
		setwd("/project/inertia_project/Data")
	} else {
		setwd("C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Data")
	}

	plans <- get(load("ca_plan_characteristics_AUG032019_small")) # plan characteristics (includes same plans in different rating areas and years)
	geog_factors <- read.csv("ca_cms_geographic_factors.csv",header=T,row.names=1)
	
# Clean-up Steps
	
	# Now we can drop the "uninsured" column from the choices and premiums object, and the uninsured row from the plan object
	plans <- plans[!is.na(plans$region),]

	# Drop CSR plans
	plans <- plans[!plans$AV %in% c(0.73,0.87,0.94),]

	# Drop 2019
	plans <- plans[plans$Year < 2019,]

	# Drop Catastrophic
	plans <- plans[plans$Metal_Level != "Minimum Coverage",]
	
	# Remove Duplicate Bronze and Small Plans
	plans$region_year <- paste(plans$region,plans$Year,sep="_")
	plans$unique_id_nohsacat <- paste(plans$Plan_Name_Small_NOHSACAT,plans$Year,plans$region,sep="_")	
	duplicated_plans <- plans[which(duplicated(plans$unique_id_nohsacat)),"unique_id_nohsacat"]
	for(mkt_yr in unique(plans$region_year)) {
		mkt_duplicated_plans <- intersect(duplicated_plans,plans[plans$region_year == mkt_yr,"unique_id_nohsacat"])
		for(j in mkt_duplicated_plans) {
			plans[plans$region_year == mkt_yr & plans$unique_id_nohsacat == j,"Premium"] <- min(plans[plans$region_year == mkt_yr & plans$unique_id_nohsacat == j,"Premium"])
			plans[plans$region_year == mkt_yr & plans$unique_id_nohsacat == j,"Hausman_instrument"] <- min(plans[plans$region_year == mkt_yr & plans$unique_id_nohsacat == j,"Hausman_instrument"])
		}
	}	
	plans <- plans[!duplicated(plans$unique_id_nohsacat),]
	
# Define covariates
	
	year_FEs <- paste("Year",c(2015:2018),sep="") 
	market_FEs <- paste("Rating_Area",c(2:19),sep="")
	
	if(add_year_FEs & add_market_FEs) {
		covariates <- c(product_covariates,year_FEs,market_FEs,instrumental_variables)
	} else if(add_year_FEs) {
		covariates <- c(product_covariates,year_FEs,instrumental_variables)
	} else if(add_year_FEs) {
		covariates <- c(product_covariates,market_FEs,instrumental_variables)
	} else {
		covariates <- c(product_covariates,instrumental_variables)
	}
	
# Create Variables
		
	# HMO Dummy (Lump EPO/PPO into one)
	plans$HMO <- as.numeric(plans$PLAN_NETWORK_TYPE == "HMO")
	plans$Anthem_HMO <- plans$Anthem * plans$HMO
	
	# Year dummies
	for(yr in c(2015:2018)) {
		plans[,paste("Year",yr,sep="") ] <- as.numeric(plans$Year == yr)
	}
		
	# Geographic Cost Factors
	plans$geog_factor <- geog_factors[paste(plans$region,plans$Year,sep="_"),"GCF"]
	
	# Rating Area Dummies
	for(mkt in c(2:19)) {
		plans[,paste("Rating_Area",mkt,sep="")] <- as.numeric(plans$Rating_Area == mkt)
	}
	
	# Rating Areas
	plans$Rating_Area <- as.factor(plans$Rating_Area)
	
	plans$Premium <- plans$Premium/1.278 # convert from 40-year-old to 21-year-old premium
	
	fields <- c("Plan_Name_Small_NOHSACAT","Year","Rating_Area","Insurer","Metal_Level",product_covariates,instrumental_variables)
	plans_to_save <- plans[,fields]
	colnames(plans_to_save)[1] <- "Plan_Name"
	rownames(plans_to_save) <- paste(plans_to_save$Plan_Name,plans_to_save$Year,plans_to_save$Rating_Area,sep="_")
	write.csv(plans_to_save,"plan_instrument_file.csv")
	
	
	
	
	
	
	
	
# Estimate Reduced form model

	
	#spec <- Premium ~ AV + Silver + HMO + Anthem + Blue_Shield + Kaiser + Health_Net + Anthem_HMO + 
	#	geog_factor + HMO_own + HMO_other + AV_own + AV_other + Hausman_instrument + as.factor(Year) + as.factor(Rating_Area)
	#reduced_form <- lm(spec,data=plans)
	reduced_form <- lm(Premium ~ .,data=plans[,covariates]) 
	summary(reduced_form)
	save(reduced_form,file="reduced_form")
	write.csv(summary(reduced_form)[4],paste("reduced_form_results.csv",sep=""))
	plans$residual <- residuals(reduced_form)
	gc()
	
	
	
#### Save Residuals

	rownames(plans) <- paste(plans$Plan_Name_Small_NOHSACAT,plans$Year,plans$region,sep="")
	if(network_flag) {
		households <- get(load("ca_household_characteristics_JUN182021_net")) # household characteristics (I * K1)
		julia_output <- read.csv(file="ca_julia_data_JUN182021_net.csv",header=TRUE)
	} else if(churn_flag) {
		households <- get(load("ca_household_characteristics_AUG032019_churn")) # household characteristics (I * K1)
		julia_output <- read.csv(file="ca_julia_data_AUG032019_churn.csv",header=TRUE)
	} else {
		households <- get(load("ca_household_characteristics_AUG032019_small")) # household characteristics (I * K1)
		julia_output <- read.csv(file="ca_julia_data_AUG032019_small.csv",header=TRUE)
	}
	julia_output$rating_area <- households[julia_output$household_number,"rating_area"]
	hrefs <- paste(julia_output[julia_output$year < 2019 & julia_output$plan_name != "Uninsured","plan_name"],
				   julia_output[julia_output$year < 2019 & julia_output$plan_name != "Uninsured","year"],
				   julia_output[julia_output$year < 2019 & julia_output$plan_name != "Uninsured","rating_area"],sep="")
	
	julia_output$residual <- 0
	julia_output[julia_output$year < 2019 & julia_output$plan_name != "Uninsured","residual"] <- plans[hrefs,"residual"]
	if(network_flag) {
		write.csv(julia_output,file="ca_julia_data_JUN182021_net.csv",row.names=FALSE)			
	} else if(churn_flag) {
		write.csv(julia_output,file="ca_julia_data_AUG032019_churn.csv",row.names=FALSE)
	} else{
		write.csv(julia_output,file="ca_julia_data_AUG032019_small.csv",row.names=FALSE)	
	}
	
	
	
	
	